function [  ] = changeDetection( varargin )
% A function for image differcing.

% basePic(the base picture) and chanPic(the warp picture) must be pre-processed.
% You can provide their path as basePicPath and chanPicPath.

thresh = -1;

if (nargin == 2)
    basePicPath = varargin{1};
    chanPicPath = varargin{2};
    detectMethod = 'Manual';
elseif (nargin == 3)
    if (isnumeric(varargin{3})),
        basePicPath = varargin{1};
        chanPicPath = varargin{2};
        detectMethod = 'Manual';
        thresh = varargin{3};
    else
        basePicPath = varargin{1};
        chanPicPath = varargin{2};
        detectMethod = varargin{3};
    end 
else
    fprintf('The number of input is wrong!\n');
    return
end

[basePic, ] = loadEnviData(basePicPath);
[chanPic, ] = loadEnviData(chanPicPath);

basePic = permute(basePic, [2, 1, 3]);
chanPic = permute(chanPic, [2, 1, 3]);

diffPic = chanPic - basePic;
[h, w, d] = size(diffPic);
diffPic = reshape(diffPic, [h * w, d]);
mag = sqrt(sum(diffPic.^2, 2));

if (strcmp(detectMethod, 'Kmeans') == 1),
    cd = kmeans(mag, 2) - 1;
    cd = reshape(cd, [h, w]);
elseif (strcmp(detectMethod, 'Otsu') == 1),
    thresh = otsuMode(mag);
    cd = (mag > thresh);
    cd = reshape(cd, [h, w]);
elseif (strcmp(detectMethod, 'KI') == 1),
    maxMag = max(mag);
    minMag = min(mag);
    mag = (mag - minMag) / (maxMag - minMag) * 255;
    mag = uint8(floor(mag));
%     mag = reshape(mag, [h, w]);
    cd = KIMode(mag);
    cd = reshape(cd, [h, w]);
elseif (strcmp(detectMethod, 'Manual') == 1),
    if (thresh == -1)
        thresh = manualMode(mag);
    end
    cd = (mag > thresh);
    cd = reshape(cd, [h, w]);
else
    fprintf('You can separate CHANGED with NON-CHANGED using as ''Kmeans'', ''Otsu'', ''KI'', or ''Manual''.\n');
    return
end

figure('NumberTitle', 'off', 'Name', ['Change Detection Using ', detectMethod, ' Method']);
imshow(cd);
end

function thresh = otsuMode(mag)
% Otsu method.
figure('NumberTitle', 'off', 'Name', 'Histogram')
hist(mag, 100);

maxMag = max(mag);
minMag = min(mag);
mag = (mag - minMag) / (maxMag - minMag);
thresh = graythresh(mag) * (maxMag - minMag) + minMag;


end

function cd = KIMode(imag)
% KITTLERMET binarizes a gray scale image 'imag' into a binary image
% Input:
%   imag: the gray scale image, with black foreground(0), and white
%   background(255).
% Output:
%   imagBW: the binary image of the gray scale image 'imag', with kittler's
%   minimum error thresholding algorithm.

% Reference:
%   J. Kittler and J. Illingworth. Minimum Error Thresholding. Pattern
%   Recognition. 1986. 19(1):41-47

MAXD = 100000;
imag = imag(:,:,1);
[counts, x] = imhist(imag);  % counts are the histogram. x is the intensity level.
GradeI = length(x);   % the resolusion of the intensity. i.e. 256 for uint8.
J_t = zeros(GradeI, 1);  % criterion function
prob = counts ./ sum(counts);  % Probability distribution
meanT = x' * prob;  % Total mean level of the picture
% Initialization
w0 = prob(1);   % Probability of the first class
miuK = 0;   % First-order cumulative moments of the histogram up to the kth level.
J_t(1) = MAXD;
n = GradeI-1;
for i = 1 : n
    w0 = w0 + prob(i+1);
    miuK = miuK + i * prob(i+1);  % first-order cumulative moment
    if (w0 < eps) || (w0 > 1-eps)
        J_t(i+1) = MAXD;    % T = i
    else
        miu1 = miuK / w0;
        miu2 = (meanT-miuK) / (1-w0);
        var1 = (((0 : i)'-miu1).^2)' * prob(1 : i+1);
        var1 = var1 / w0;  % variance
        var2 = (((i+1 : n)'-miu2).^2)' * prob(i+2 : n+1);
        var2 = var2 / (1-w0);
        if var1 > eps && var2 > eps   % in case of var1=0 or var2 =0
            J_t(i+1) = 1+w0 * log(var1)+(1-w0) * log(var2)-2*w0*log(w0)-2*(1-w0)*log(1-w0);
        else
            J_t(i+1) = MAXD;
        end
    end
end
minJ = min(J_t);
index = find(J_t == minJ);
th = mean(index);
th = (th-1)/n
cd = im2bw(imag, th);

% figure, imshow(imagBW), title('kittler binary');
end

function thresh = manualMode(mag)
% Manual choose.
figure('NumberTitle', 'off', 'Name', 'Histogram')
hist(mag, 100);
title('Click in the histogram to set a threshold.');

maxMag = max(mag);
minMag = min(mag);
thresh = -1;
while (thresh > maxMag || thresh < minMag),
    [thresh, ~] = ginput(1);
end
fprintf('You chose %g as the threshold.\n', thresh);
% mag = reshape(mag, [h, w]);
% cd = (mag > thresh);
% figure(2)
% imshow(cd);

end
